** Siyaset Biliminde Etkileşimli Modeller Üzerine: Yineleme Kodu 12.02.22
cd "/Users/mmoral/Dropbox/Sabancı University/-ONGOING/Siyaset Biliminde Etkileşimli Modeller Üzerine/"

clear

set seed 110222							/* Rastlantısal olarak türetilen değişkenler ve bunlar kullanılarak tahmin edilen regresyonların makalede raporlanan değerleri/katsayıları vermesi için aynı `seed' numarası kullanılmalıdır. */ 

/*
foreach prog in ftest grc1leg2 interflex estout mplotoffset moremata{
cap net install `prog'.pkg, replace		/* Grafik ve hesaplamalarda kullanılan ilave programlar... */
cap ssc install `prog', replace			/* ...ayrıca yüklenmedikleri takdirde yineleme kodunun bir bölümü çalışmayacaktır. */
}
*/

set scheme mmoral5						/* Grafiklerde kullanılan şema yazardan temin edilebilir. */

** ETKİLEŞİMLİ DOĞRUSAL MODELLER
set obs 150								/* Şekil 1'de iki farklı grup gözlemin ayırt edilebilmesi amacıyla gözlem sayısı (N), 150 olarak belirlenmiştir. */

gen x=2*runiform()						/* 0-2 arasında değerler alacak, tekdüze dağılan rastlantısal X değişkeni türetilmektedir. */
gen z=(runiform()>=.5)					/* Yine rastlantısal olarak türetilen ve tekdüze dağılan sürekli değişkenin değeri 0.5'e eşit ya da 0.5'ten büyükse, iki değerli Z değişkeni 1 değerini almaktadır. */
gen e=rnormal()							/* Sentetik veriler türetilirken hataların standart normal dağıldığı varsayımı yapılmıştır. */

gen y=1+z+x*z+e 						/* Y değişkeni bu eşitlik ile türetilmektedir. */

* Toplamlı ve Etkileşimli Doğrusal Modellerin Karşılaştırması
eststo m1a: reg y x z
*reg y x z c.x#c.z 						/* Etkileşimli doğrusal model, Stata'nın operatörleri kullanılarak da tahmin edilebilir. "#" sembolü çarpımsal etkileşim, "c." makrosu da sürekli değişkenler için kullanılmaktadır. */
eststo m1i: reg y c.x##c.z 				/* "##" operatörü kullanıldığındaysa kurucu değişkenlerin ayrı ayrı belirtilmesine gerek kalmamaktadır. */

ftest m1a m1i							/* F-test yardımıyla toplamlı ve etkileşimli modeller karşılaştırılarak, etkileşimli modelin veriye uyumunun aynı ya da daha düşük olduğuna dair boş hipotez test edilebilir. */

** Şekil1: Toplamlı ve Etkileşimli Doğrusal Modeller
est restore m1a 						/* Tahmin edilen Y değerleri, yukarıda m1a olarak adlandırılan toplamlı model kullanılarak hesaplanmaktadır. */

estat sum								/* Bağımlı ve bağımsız değişkenlerin örneklem içerisinde aldıkları değerler özetlenmektedir. */

margins, at(x=(0(.1)2) z=(0 1))			/* Eşitlikteki her iki kurucu değişken için margins komutu aracılığıyla Y'yi tahmin etmekte kullanılacak değerleri belirlenmektedir. */

marginsplot, recastci(rarea) recast(line) ci1opts(fcolor(gs12%65) lwidth(none)) ci2opts(fcolor(gs6%65) lwidth(none)) ytitle(Y) xtitle(X) title("Şekil 1a. Toplamlı Doğrusal Model (1T)") ///
addplot((scatter y x if z==0 & e(sample), msymbol(oh) msize(small) mlwidth(thin) mcolor(black%25)) ///
(scatter y x if z==1 & e(sample), msymbol(plus) msize(small) mlwidth(thin) mcolor(black%75) xsca(r(0 2.5)) xlab(0(.5)2))) ///
text (3.58 2.25 "`=ustrunescape("Y\u0302")'={&beta}{subscript:0}+{&beta}{subscript:2}+({&beta}{subscript:1})X", just(left) size(vsmall)) ///
text (1.51 2.2 "`=ustrunescape("Y\u0302")'={&beta}{subscript:0}+({&beta}{subscript:1})X", just(left) size(vsmall)) name(Fig1a, replace)

est restore m1i 						/* Tahmin edilen Y değerleri, yukarıda m1i olarak adlandırılan etkileşimli model kullanılarak hesaplanmaktadır. */

margins, at(x=(0(.1)2) z=(0 1))			/* Aynı örneklem kullanıldığından, Y'nin değerinin tahmin edilmesinde X ve Z'nin aynı değerleri kullanılmaktadır. */

marginsplot, recastci(rarea) recast(line) ysca(alt) ci1opts(fcolor(gs12%65) lwidth(none)) ci2opts(fcolor(gs6%65) lwidth(none)) ytitle(Y) xtitle(X) title("Şekil 1b. Etkileşimli Doğrusal Model (1E)") ///
addplot((scatter y x if z==0 & e(sample), msymbol(oh) msize(small) mlwidth(thin) mcolor(black%25)) (scatter y x if z==1 & e(sample), msymbol(plus) msize(small) mlwidth(thin) mcolor(black%75) xsca(r(0 2.5)) xlab(0(.5)2))) ///
text (4.3 2.3 "`=ustrunescape("Y\u0302")'={&beta}{subscript:0}+{&beta}{subscript:2}+({&beta}{subscript:1}+{&beta}{subscript:3})X", just(left) size(vsmall)) ///
text (.85 2.175 "`=ustrunescape("Y\u0302")'={&beta}{subscript:0}+({&beta}{subscript:1})X", just(left) size(vsmall)) name(Fig1b, replace)

grc1leg2 Fig1a Fig1b, ycommon legendfrom(Fig1b) 					/* Bu iki grafik, y aksı ortak olacak şekilde birleştirilmektedir. */

gr_edit .legend.plotregion1.label[1].text={}						/* Bu bölümde Şekil 1'in lejandları düzenlenmektedir. */
gr_edit .legend.plotregion1.label[1].text.Arrpush Ŷ | Z=0
gr_edit .legend.plotregion1.label[2].text={}
gr_edit .legend.plotregion1.label[2].text.Arrpush Ŷ | Z=1
gr_edit .legend.plotregion1.label[3].text={}
gr_edit .legend.plotregion1.label[3].text.Arrpush Y | Z=0
gr_edit .legend.plotregion1.label[4].text={}
gr_edit .legend.plotregion1.label[4].text.Arrpush Y | Z=1

gr export "Şekil1.jpg", as(jpg) quality(100) replace

** Çarpımlı Etkileşimli ve Toplamlı Doğrusal Modeller
eststo m2: reg y x if z==0 				/* Z'nin gözlemlenmediği alt örneklemde X'in Y'ye doğrudan etkisi tahmin edilmektedir. */
eststo m3: reg y x if z==1 				/* Z'nin gözlemlendiği alt örneklemde X'in Y'ye doğrudan etkisi tahmin edilmektedir. */

** Şekil2: X ve Z'nin Y'ye Marjinal Etkileri (Etkileşimli Doğrusal Model)
est restore m1i

margins, dydx(x) at(z=(0 1)) 			/* Margins komutu yardımıyla Y'nin X'e göre birinci mertebeden türevi alınmakta ve X'in Y'ye marjinal etkisi Z'nin her iki değeri için hesaplanmaktadır. */
mat b_est=r(b)

marginsplot, recast(scatter) recastci(rspike) xsca(r(-.125 1.125)) ///
text (-0.22 0.0375 "{&beta}{subscript:1}", just(left) size(vsmall)) ///
text (1.17 .925 "{&beta}{subscript:1}+{&beta}{subscript:3}", just(left) size(vsmall)) ///
title("Şekil 2a. X'in Y'ye Marjinal Etkisi | Z") xtitle(Z) ytitle("Marjinal Etki") yline(0, lcolor(red) lwidth(thin)) name(marjx, replace)

gen nerede=-1							/* Halı çizişinin grafikte gösterileceği Y değeri belirlenmektedir. */
gen str tik="|"							/* Halı çizişinde kullanılacak sembol, nümerik olmayan bir değişken olarak türetilmektedir. */

margins, dydx(z) at(x=(0(.05)2)) 
marginsplot, recastci(rspike) plotopts(mlwidth(none)) ciopts(lwidth(thin)) xsca(r(0 2.25)) ///
title("Şekil 2b. Z'in Y'ye Marjinal Etkisi | X") xtitle(X) ytitle("Marjinal Etki") ///
text (3.45 2.125 "{&beta}{subscript:2}+{&beta}{subscript:3}X", just(left) size(vsmall)) ///
addplot((scatter nerede x if x, ms(i) mlab(tik) mlabsize(vsmall) mlabgap(-1) mlabpos(6))) yline(0, lcolor(red) lwidth(thin)) legend(off) name(marjz, replace)

grc1leg2 marjx marjz, loff legendfrom(marjx) ytol1title ycommon		/* Bu iki grafik, y aksının başlığı ortak olarak şekilde birleştirilmektedir. */

gr export "Şekil2.jpg", as(jpg) quality(100) replace

* Ortalama Marjinal Etkilerin Hesaplanması
tab z
disp (b_est[1,1]*71+b_est[1,2]*79)/`r(N)'							/* X'in ortalama marjinal etkisi bu şekilde hesaplanmaktadır */
margins, dydx(x)													/* Aynı değer, dydx opsiyonu ile de hesaplanabilecektir. */

* Marjinal Etkilerin Manuel Olarak Hesaplanması (Dipnot#7 ve Dipnot#8)
mat b=e(b)								/* Margins komutuna alternatif olarak, marjinal etkiler beta ve varyans-kovaryans matrislerinden manuel olarak da hesaplanabilir. */
mat V=e(V)								/* Bu amaçla ilgili matrislerden gerekli sayıllar elde edilmelidir. */
sca b2=b[1,2]							/* Katsayılar */
sca b3=b[1,3]
sca varb2=V[2,2]						/* Varyanslar */
sca varb3=V[3,3]
sca covb2b3=V[2,3]						/* Kovaryans */

gen fakex=(_n-1)*0.05 in 1/41			/* 0 ile 2 arasında 0.05 birim kadar artan değerler alan sahte bir değişken türetilerek (2/.05+1 sayıda gözlem), ilgili marjinal etkiler ve standart hataları hesaplanmaktadır. */ 
gen marg_zhat=b2+b3*fakex											/* Marjinal etki */	
gen se_zhat=sqrt(varb2+fakex^2*varb3+2*fakex*covb2b3)				/* Standart hata */
gen ub_zhat=1.976*se_zhat+marg_zhat									/* Grafikte kullanılacak güven aralığı oluşturulmaktadır. */
gen lb_zhat=-1.976*se_zhat+marg_zhat 								/* Bu amaçla kullanılacak kritik değer, hataların serbestlik derecesine ve tercih edilen güven aralığına göre t tablosundan seçilmelidir. */
gen t_zhat=marg_zhat/se_zhat

list fakex *zhat* in 1/41				/* Manuel ve otomatik (margins komutuyla) olarak hesaplandığında aynı marjinal etki tahminlerinin elde edildiği görülecektir. */

** Şekil3: X'in Y_2'ye Doğrusal Olmayan Marjinal Etkisi | Z_2 (Etkileşimli Doğrusal Model
gen z2=rnormal()						/* Standart normal dağılan, sürekli Z2 değişkeni türetilmektedir. */

gen y2=0.5 -1.5*x*z2 + 2.5*x*z2^2 + e	/* Y2 değişkeni, bu eşitlikten türetilmektedir */

eststo m4: reg y2 c.x##c.z2

egen z2_gr=cut(z2), group(3)			/* Z2 değişkeni, interflex yönteminin açıklanması ve grafikte yer alan toplamlı modeller için eşit sayıda gözlemden oluşan, eşit sayıda gözleme sahip üç gruba bölünmektedir. */

twoway (scatter y2 x, msize(small) msymbol(+) mlwidth(thin) mcolor(black%75)) ///
(lfitci y2 x if z2_gr==0, fcolor(maroon%75) alwidth(none) lcolor(white%50)) ///
(lfitci y2 x if z2_gr==1, fcolor(dkorange%75) alwidth(none) lcolor(white%50)) ///
(lfitci y2 x if z2_gr==2, fcolor(dkgreen%75) alwidth(none) lcolor(white%50)), ///
xtitle(X) ytitle(Y{subscript:2}) title("Şekil 3a. Y{subscript:2}'nin Gözlemlenen ve Tahmin Edilen Değerleri") ///
legend(ring(0) rows(1) size(vsmall) pos(11) order(2 "Ŷ{subscript:2} | Z{subscript:2} - D(üşük)" 4 "Ŷ{subscript:2} | Z{subscript:2} - O(rta)" 6 "Ŷ{subscript:2} | Z{subscript:2} - Y(üksek)")) name(Şekil3a, replace)

* Interflex ve Grafiklenecek İstatistiklerin Elde Edilmesi
interflex y2 x z2, title("Şekil 3b. Interflex") xlabel(Z{subscript:2}) dlabel(X) ylabel(Y{subscript:2}) bw(.1)	/* Interflex aracının görselleştirdiği histogramların bant aralıkları değiştirilebilmektedir. */

mat me=r(margeff)						/* Margins komutunda kullanılan değişken değerleri `r(at)' matrisinden, tahmin edilen değerleri içeren tablo ise `r(table)' matrisinden elde edilebilir. */ 
mat bin=r(estBin)						/* `r(table)' matrisi ilgili değişkenleri türetebilmek için transpoze edilmektedir. */
svmat me, name(etki)					/* Her iki matristeki kolonlardan ayrı değişkenler üretilmektedir. */
svmat bin, name(bin)
gen str bin6=""
gen bin7=.
tokenize "D O Y"
forval i=1/3{
replace bin6="``i''" in `i'
replace bin7=9.25 in `i'
}

twoway (line etki2 etki1, lwidth(thin) lcolor(black%75)) (rarea etki4 etki5 etki1, lwidth(none) fcolor(gs12%65)) ///
(rcap bin4 bin5 bin1, lcolor(maroon%50) lwidth(thin)) (scatter bin2 bin1, mcolor(maroon%75) msymbol(o) mlwidth(none)) ///
(scatter bin7 bin1, mcolor(gs12%65) msymbol(i) mlab(bin6) mlabpos(12) mlabsize(medlarge)) (hist z2 if e(sample), width(.1) lwidth(vthin) fcolor(gs12%65) yaxis(2) percent) ///
, yline(0, lwidth(vthin) lcolor(red)) xtitle("Z{subscript:2}") ytitle("X'in Y{subscript:2}'ye Marjinal Etkisi") title("Şekil 3b. Interflex") ylab(0(4)28, axis(2))  ysca(axis(2) off) legend(off) name(Şekil3b, replace)

gr combine Şekil3a Şekil3b				/* Bu iki grafik birleştirilmektedir. */

gr export "Şekil3.jpg", as(jpg) quality(100) replace

* Interflex Yönteminin Detaylı Açıklaması 
interflex y2 x z2, title("Şekil 3b. Interflex") xlabel(Z{subscript:2}) dlabel(X) ylabel(Y{subscript:2}) bw(.1)	/* Interflex aracının görselleştirdiği histogramların bant aralıkları değiştirilebilmektedir. */

mat lis r(estBin)						/* Interflex'in kullandığı üç alt örneklem, bu örneklemlerdeki ortanca değerler ve bu değerler için hesaplanan marjinal etkiler komut çalıştırıldıktan sonra ilgili matristen elde edilebilir. */ 
tabstat z2, by(z2_gr) stats(p50)		/* Bu değerler ile yukarıda oluşturulan grupların ortanca değerlerinin aynı olduğu görülecektir. */

loc nbin=3								/* İlk olarak kaç alt örneklem oluşturulacağı belirlenmelidir. */
egen cutz=cut(z2), gr(`nbin')			/* Ardından koşullandırıcı değişken bu sayıda alt gruba bölünmekte, */
tab cutz, gen(gr)						/* bundan bir az sayıda grup için kukla değişkenler oluşturulmakta, */
forval i=2/`nbin'{						/* ilk alt örneklem dışındaki diğer alt örneklemler için */
qui sum z2 if gr`i'==1, detail			/* bu alt grupların ortanca noktaları belirlenerek, */	
gen int`i'_1=gr`i'*x 					/* kukla değişkenler ile her iki kurucu değişken ve bunların çarpımlı etkileşimlerin çarpımlı etkileşimleri oluşturulmaktadır. */
gen int`i'_2=gr`i'*(z2-`r(p50)') 		/* Z2'nin bu iki alt örneklemde aldığı değerler, ortanca noktalar 0'a eşit olacak şekilde yeniden ölçeklendirilerek belirlenmelidir. */
gen int`i'_3=gr`i'*x*(z2-`r(p50)')
}

eststo mm1: reg y2 c.x##c.z2 int2_1 gr2 int2_2 int2_3 int3_1 gr3 int3_2 int3_3			/* Tüm bu terimleri içeren doğrusal model ile interflex yöntemi kullanılarak tahmin edilen doğrusal model aynı katsayıları verecektir. . */
eststo mm2: interflex y2 x z2

return list
mat lis r(estBin)						/* Interflex yöntemi ile tahmin edilen ve Şekil 3b'de gösterilen marjinal etki (nokta) tahminleri ile bunların Z2'nin hangi değerleri için hesaplandığı `estbin' isimli matriste saklanmaktadır. */

esttab mm1 mm2, se stats(N r2)			/* İki model karşılaştırıldığında, katsayıların nasıl hesaplandığı görülecektir. */
 
reg y2 c.x##c.z2 if gr2==0 & gr3==0		/* Etkileşimli model ilk olarak birinci alt örneklemimiz için (Z2'in düşük değerleri) için tahmin edilmektedir. */
margins, dydx(x) at(z2=(-.87722316))	/* X'in bu alt örneklemdeki marjinal etkisi, Y2'nin X'e göre birinci mertebeden türevi alınarak, Z2'nin bu alt örneklemdeki ortanca noktası için hesaplanmaktadır. . */
reg y2 int2_1 int2_2 int2_3	if gr2		/* Diğer alt örneklemler için (Z2 değişkeninden Z2'nin bu alt örneklemlerdeki ortanca değeri çıkarıldığından) X'in bir birimlik değişimi (0-1) Z2'nin ortanca noktasındaki marjinal etkisini göstermektedir. */
disp e(b)[1,1]
reg y2 int3_1 int3_2 int3_3	if gr3
disp e(b)[1,1]

** Tablo1 - Toplamlı ve Etkileşimli Doğrusal Modeller
esttab m1a m1i m2 m3 m4 using Tablo1.csv, replace b(%10.3f) se stats(N r2, fmt(0 3) labels("N" "\$R^{2}\$")) ///
unstack starlevels(* 0.1 ** 0.05 *** 0.01) alignment(l) nonote nogaps nonumbers noomit nobase compress mtitle("Model.1T" "Model.1E" "Model.2" "Model.3" "Model.4") ///
coeflabels(_cons Constant) substitute(\_ _) fragment

** ETKİLEŞİMLİ DOĞRUSAL-OLMAYAN MODELLER
clear				
set obs 1000							/* Etkileşimli doğrusal-olmayan model örnekleri için 1000 gözlemli yeni bir veriseti oluşturulmaktadır. */

gen x1=rnormal()						/* Standart normal dağılan x1, tekdüze dağılan x2 ve iki-değerli x3 değişkenleri oluşturulmaktadır. */
gen x2=runiform()
gen x3=rbinomial(1,.3)

gen xb=0.25+.5*x1*x2+ x1*x2^2+ .5*x3+rnormal()											/* Bu değişkenler ve x1 ile x2'nin çarpımlı etkileşimleri ve standart normal dağılan hatalar kullanılarak sistematik unsur oluşturulmakta, */  
gen pr=1/(1+exp(-xb))																	/* ardından lojistik kümülatif dağıtım fonksiyonu ile bu değerler olasılık uzayına aktarılmakta, */
gen y3=pr<.5 																			/* son olarak ilgili gözlemlenme olasılıklarının .5'den yüksek olduğu değerler için 1 değerini alan iki-değerli Y3 değişkeni üretilmektedir. */

* Olabilirlik Oranı Testi
eststo mm3: logit y3 x1 x2 x3															/* x1, x2 ve x3'ün bağımsız değişkenler olduğu etkileşimli lojistik regresyon tahmin edilmektedir. */
eststo mm4: logit y3 c.x1##c.x2 x3														/* x1, x2, (x1 x x2) ve x3'ün bağımsız değişkenler olduğu etkileşimli lojistik regresyon tahmin edilmektedir. */
lrtest mm4 mm3

* Margins Komutu ile Y3*'ün Beklenen Değerlerinin Hesaplanması
sum x1, detail							/* x1 değişkeni özetlenmektedir. */
margins, at(x1=(`r(min)' `r(mean)' `=`r(mean)'+1' `r(max)') x2=(0(.1)1) x3=0)			/* x1 değişkeninin asgari, ortalama, ortalama+1 birim ve azami değerlerinin her biri için x2 değişkeni örneklem içerisindeki aralığı içerisindeki (0-1) her nokta olasılık tahmini için ilgili değer 0.1 birim arttırılmakta, x3 değişkeni ise modunda/ortanca değerinde sabit tutularak Y3*'ün gözlemlenme olasılığı (Pr(Y3=1)) hesaplanmaktadır. */	

* Marginsplot Komutu
marginsplot, xdim(x2) recastci(rspike) legend(order(5 "X{subscript:1}=R{subscript:min.}" 6 "X{subscript:1}=R{subscript:ort.}" 7 "X{subscript:1}=R{subscript:ort.+1}" 8 "X{subscript:1}=R{subscript:maks.}")) ///
ytitle(Pr(Y{subscript:3}=1)) xtitle(X{subscript:2}) title("Şekil 4a. Y{subscript:3}'nin Gözlemlenme Olasılığı") 	/* Hesaplanan olasılıklar marginsplot komutu ile de görselleştirilebilir. */

* Mplotoffset Paketi
mplotoffset, offset(.008) xdim(x2) recastci(rspike) ciopts(lwidth(thin)) plotopts(msymbol(o) mlwidth(none) msize(small) lpattern(shortdash) lwidth(thin)) ///
ci1opts(lcolor(black%50)) ci2opts(lcolor(maroon%50)) ci3opts(lcolor(dkorange%50)) ci4opts(lcolor(dkgreen%50)) ///
plot1opts(lcolor(black%75) mcolor(black%75)) plot2opts(lcolor(maroon%75) mcolor(maroon%75)) plot3opts(lcolor(dkorange%75) mcolor(dkorange%75)) plot4opts(lcolor(dkgreen%75) mcolor(dkgreen%75)) ///
legend(order(5 "X{subscript:1}=R{subscript:min.}" 6 "X{subscript:1}=R{subscript:ort.}" 7 "X{subscript:1}=R{subscript:ort.+1}" 8 "X{subscript:1}=R{subscript:maks.}")) ///
ytitle(Pr(Y{subscript:3}=1)) xtitle(X{subscript:2}) title("Şekil 4a. Y{subscript:3}'ün Gözlemlenme Olasılığı") name(Şekil4a, replace)

* Margins Komutundan İlgili İstatistiklerin Elde Edilmesi
mat at=r(at)							/* Margins komutunda kullanılan değişken değerleri `r(at)' matrisinden, tahmin edilen değerleri içeren tablo ise `r(table)' matrisinden elde edilebilir. */ 
mat r=r(table)'							/* `r(table)' matrisi ilgili değişkenleri türetebilmek için transpoze edilmektedir. */
svmat at, name(ats)						/* Her iki matristeki kolonlardan ayrı değişkenler üretilmektedir. */
svmat r, name(pred)

egen gr=group(ats1) if !missing(ats1)	/* x1 değişkeni için margins komutunda atanan dört değere tekabül eden dört ayrı grup oluşturulmaktadır. */
gen ats2_v2=ats2-.024+gr*.008			/* Güven aralıklarının grafikte ayırt edilebilmesi için jitter eklenmektedir. */

* Marjinal Etki Tahminlerinin Manuel Olarak Görselleştirilmesi
twoway (rspike pred5 pred6 ats2_v2 if gr==1, lcolor(black%50) lwidth(thin)) (connected pred1 ats2_v2 if gr==1, msymbol(o) mlwidth(none) msize(small) color(black%75) lpattern(shortdash) lwidth(thin)) ///
(rspike pred5 pred6 ats2_v2 if gr==2, lcolor(maroon%50) lwidth(thin)) (connected pred1 ats2_v2 if gr==2, msymbol(o) mlwidth(none) msize(small) color(maroon%75) lpattern(shortdash) lwidth(thin)) ///
(rspike pred5 pred6 ats2_v2 if gr==3, lcolor(dkorange%50) lwidth(thin)) (connected pred1 ats2_v2 if gr==3, msymbol(o) mlwidth(none) msize(small) color(dkorange%75) lpattern(shortdash) lwidth(thin)) ///
(rspike pred5 pred6 ats2_v2 if gr==4, lcolor(dkgreen%50) lwidth(thin)) (connected pred1 ats2_v2 if gr==4, msymbol(o) mlwidth(none) msize(small) color(dkgreen%75) lpattern(shortdash) lwidth(thin)), ///
legend(order(2 "X{subscript:1}=R{subscript:min.}" 4 "X{subscript:1}=R{subscript:ort.}" 6 "X{subscript:1}=R{subscript:ort.+1}" 8 "X{subscript:1}=R{subscript:maks.}")) ///
ytitle(Pr(Y{subscript:3}=1)) xtitle(X{subscript:2}) title("Şekil 4a. Y{subscript:3}'ün Gözlemlenme Olasılığı")

* Simülasyon Yöntemiyle Bir Sürekli Değişkenin (Ortalama Değerindeki) Marjinal Etkisinin Hesaplanması
preserve
logit
estat sum
qui sum x1
loc x1_ort=`r(mean)'																	/* İlk farkı alınacak değerlerin hesaplamasında x1 değişkeninin örneklem içerisindeki ortalama değeri kullanılmaktadır. */

drawnorm b1-b5, n(10000) means(e(b)) cov(e(V)) clear									/* beta ve varyans-kovaryans matrisleri kullanılarak 10,000 set (çok terimli normal dağılan) regresyon katsayısı simüle edilmektedir. */

postutil clear
postfile mm ola0 as0 us0 ola1 as1 us1 olafark asfark usfark using sim, replace			/* Basit bir program yardımıyla, x2'nin aldığı her farklı değer için bu satırda listelenen tüm değişkenler `sim' ismindeki verisetine kaydedilmemektedir. */ 
forval i=0/10{																			/* x1'deki bir birimlik değişimin Y3*'e marjinal etkisi x2'nin 11 farklı değeri için (simüle edilen katsayılar kullanılarak) hesaplanmaktadır. */ 
{
gen xb0=b1*(`x1_ort') 	+ b2*(`i'/10) + b3*(`x1_ort')*(`i'/10)   + b4*0 + b5			/* Her gözlem (simülasyon) için xb0 ve xb1'in (sistematik unsurlar) aldıkları değerler hesaplanmaktadır. */
gen xb1=b1*(`x1_ort'+1) + b2*(`i'/10) + b3*(`x1_ort'+1)*(`i'/10) + b4*0 + b5			
forval j=0/1{																			
gen ola`j'=1/(1+exp(-xb`j'))															/* Doğrusal sistematik unsurlar, lojistik kümülatif dağılım fonksiyonu ile olasılık uzayına aktarılmaktadır. */
egen ortolas`j'=mean(ola`j')																	
}
gen olafark=ola1-ola0																	/* xb1 ve xb0 arasındaki ilk fark hesaplanmaktadır. */
egen ortolasfark=mean(olafark)															/* Ortalama marjinal etkiler hesaplanmaktadır. */

tempname ola0 as0 us0 ola1 as1 us1 olafark asfark usfark
foreach k in 0 1 fark{
_pctile ola`k', p(2.5,97.5)																/* Güven aralıklarının alt ve üst sınırları, %95 güven aralığında belirlenmektedir. */
sca `as`k''=r(r1)
sca `us`k''=r(r2)
sca `ola`k''=ortolas`k'																	/* Simülasyonlardan elde edilen değerlerin ortalaması ise ortalama marjinal etkiyi vermektedir. */
}
post mm (`ola0') (`as0') (`us0') (`ola1') (`as1') (`us1') (`olafark') (`asfark') (`usfark')	
}
drop xb0 xb1 ola0 ola1 ortolas0 ortolas1 olafark ortolasfark
local i=`i'+1
}
postclose mm
restore

append using sim, gen(appended)										/* "sim.dta" olarak kaydedilen marjinal etkileri içeren veriseti orijinal verisetine eklenmektedir. */ 
gen sahte_x2=(_n-1001)*.1 in 1001/1011								/* Bu verisetindeki ilk satırdaki değer x2'nin asgari değerini, sonuncu ise azami değerini göstermektedir. */

twoway (hist x2 if e(sample), lwidth(vthin) fcolor(gs12%65) yaxis(2)) ///
(rspike usfark asfark sahte_x2, lcolor(gs0%95) lwidth(thin)) (connected olafark sahte_x2, lpattern(solid) msymbol(smcircle) lwidth(thin) mlwidth(none) lcolor(gs0%95) mcolor(gs0%95)), ///
legend(order(2 "Pr(Y{subscript:3}=1 | X{subscript:1}=D{subscript:ort.+1}) - Pr(Y{subscript:3}=1 | X{subscript:1}=D{subscript:ort.})"))  ///
xlab(0(.2)1) yline(0, lcolor(red) lwidth(vthin)) xtitle("X{subscript:2}") ysca(off axis(2)) ytitle("X{subscript:1}'in Ortalama Değerindeki Marjinal Etkisi") ///
title("Şekil 4b. X{subscript:1}'in Y{subscript:3}'ün Gözlemlenme Olasılığına Marjinal Etkisi") name(Şekil4b, replace)

gr combine Şekil4a Şekil4b, ycommon
gr export "Şekil4.jpg", as(jpg) quality(100) replace

* Margins Komutu Kullanılarak (Türev Yöntemiyle) Kısıtlı Bir Değişkenin (Ortalama) Marjinal Etkisinin Hesaplanması
margins, at((mean) x1 x2 x3=(0 1)) noatlegend						/* İki-değerli x3 değişkeninin her iki değeri için, süreki x1 ve x2 değişkenleri ortalama değerlerinde sabit tutularken Y3'ün gözlemlenme olasılığı hesaplanmaktadır. . */
margins, dydx(x3) at((mean) x1 x2) 									/* Bu iki olasılık arasındaki fark ile margins komutu kullanılarak kısıtlı x3'e göre türev alınması aynı (X3'ün bir birimlik değişiminin marjinal etkisi) marjinal etki tahminini verecektir. */

* Manuel Olarak ve Margins Komutu Kullanılarak (Türev Yöntemiyle) Sürekli Bir Değişkenin (Ortalama) Marjinal Etkisinin Hesaplanması
preserve
qui sum x1 if e(sample)
sca ch=(abs(r(mean))+0.0001)*0.0001									/* Stata'nın margins komutu sürekli değişkenlerin marjinal etkilerini hesaplarken ortalama değeri kullanmakta ve oldukça küçük bir değişimin ortalama etkisini hesaplamaktadır. */
clonevar x1_c=x1
qui sum x2 if e(sample)
replace x2=`r(mean)'
qui sum x3 if e(sample), detail
replace x3=`r(p50)'
replace x1=x1_c-scalar(ch)
predict lw_0 if e(sample)
replace x1=x1_c+scalar(ch)
predict lw_1 if e(sample)
gen dydx=(lw_1-lw_0)/(2*scalar(ch))
sum dydx
restore

margins, dydx(x1) at((mean) x2 (median) x3) 						/* Süreki x2 değişkeni ortalama değerinde, iki-değerli x3 ise mod/ortanca değerinde sabit tutularak, (analitik) türev alınması ile yukarıda hesapladığımız değere yakınlaşılabilecektir. */

* Margins Komutu Kullanılarak (İlk Fark Yöntemiyle) Sürekli bir Değişkenin (Ortalama Değerindeki) Marjinal Etkisinin Hesaplanması
qui sum x1 if e(sample)												/* Bu yöntemde ilk olarak ortalama değer hesaplanmakta, ardından bu değer ile bu değere 1 eklenerek hesaplanan iki gözlemlenme olasılığının farkı alınmaktadır */
margins, at((mean) x2 (median) x3 x1=(`r(mean)' `=`r(mean)'+1')) contrast(atcontrast(r._at))

drop if appended													/* Simülasyonlardan elde edilen (ortalama) marjinal etkiler, */
drop ats1-sahte_x2													/* ve yukarıdaki grafiklerde kullanılan sahte değişkenler verisetinden silinmektedir. */
erase "sim.dta"														/* Simülasyon sonuçları klasörden silinmektedir. */
tempfile DEML 														/* `DEML' isminde geçici bir veriseti oluşturularak aşağıdaki işlemlerde de aynı veriseti kullanılmaya devam etmektedir. */
save `DEML'

* Zhirnov vd. 2022
mata 																/* Yazarların kullandığı her üç program da bu bölümde tanımlanmaktadır (bkz. Zhirnov vd. (2022)). */
real matrix me(coef, x, x_new) {
dydx=logistic(coef*x_new') - logistic(coef*x')
return(dydx)
}
end

mata
real matrix me_byrow(coef, X, Z) {
dydx=me(coef, X, Z)
means=mean(dydx)'
ra=mm_quantile(dydx, 1, (0.025 \ 0.975))'
return((means,ra))
}
end

mata
real matrix me_wt(coef, X, Z, group_id, weight) {
dydx=me(coef, X, Z)
groups=uniqrows(group_id)
wtm=J(cols(dydx), rows(groups), .)
obs=J(rows(groups),1,.)
for (i=1; i<=rows(groups); i++) {
	wtmc=(group_id:==groups[i]):*weight
	obs[i]=sum(wtmc)
	wtm[.,i]=wtmc/sum(wtmc)
	}
dydxw=dydx*wtm
means=mean(dydxw)'
ra=mm_quantile(dydxw, 1, (0.025 \ 0.975))'
return((groups,obs,means,ra))
}
end

* Marjinal Etkilerin Hesaplanması
use `DEML', clear
logit y3 c.x1##c.x2 x3												/* x1, x2, (x1 x x2) ve x3'ün bağımsız değişkenler olduğu etkileşimli lojistik regresyon tahmin edilmektedir. */

mat beta=e(b)[.,e(depvar) + ":"]									/* Beta matrisinden katsayılar, varyans-kovaryans matrisinden de bunların varyans ve katsayıları elde edilmektedir. */
mat vcov=e(V)[e(depvar) + ":",e(depvar) + ":"]

preserve
drawnorm coef1-coef`=colsof(beta)', n(10000) means(beta) cov(vcov) clear 					/* Simülasyon yönteminde olduğu gibi regresyon katsayıları 10,000 kez simüle edilmekte ve bu değerler coef ismindeki matrise aktarılmaktadır. */
putmata coef=(*), replace
restore

replace x1=round(x1,0.05) 											/* 0.05 birimlik aralıklar oluşturularak, sürekli x1 ve x2'nin değerleri bu aralıklara yuvarlanmakta, bu sayede grafikte 1000 adet nokta tahmini yapılmasının önüne geçilerek, hesaplamaların daha hızlı yapılması sağlanmaktadır. Tercih edilmediği takdirde, bu bölüm çıkarılabilir. */	
replace x2=round(x2,0.05)
replace x3=0 														/* x3'ün değeri ortanca noktasında (modunda) sabit tutulmaktadır. */
collapse (mean) x3 (count) obs=x3, by(x1 x2) 						/* Grafikteki işaretçilerin büyüklükleri ortak dağılımın her bir noktasındaki gözlem sayısına orantılı olarak belirleneceğinden, her bir nokta tahmininde kaç adet gözlem olduğu sayılmaktadır. */

gen x1x2=x1*x2 														/* Çarpımlı etkileşim oluşturulmakta ve ardından ağırlıklar ile birlikte X matrisine aktarılmaktadır. */
putmata wt=obs X=(x1 x2 x1x2 x3 1), replace

preserve 
replace x1=x1+1 													/* x1'in değerine 1 eklenerek, X1 isminde yeni bir matrix oluşturulmaktadır. */
replace x1x2=x1*x2
putmata wt=obs X1=(x1 x2 x1x2 x3 1), replace
restore

mata: mebr=me_byrow(coef, X, X1)									/* Marjinal etki tahminleri, bu iki matris arasındaki ilk fark alınarak ve yukarıda tanımlanan program yardımıyla hesaplanmaktadır. */
getmata (me_est lb ub)=mebr 										/* Bu tahminlerin çevrelerindeki güven aralıklarının alt ve üst sınırları da aynı programdan elde edilmektedir. */

gen significant=(lb>0 & ub>0)|(lb<0 & ub<0)							/* İstatistiki olarak anlamlı nokta marjinal etki tahminleri işaretlenmektedir. */

gen counter=_n 														/* Grafikte kullanılacak tikler belirlenmektedir. */
qui sum counter
loc extra1=`=r(max)'+1
loc extra2=`=r(max)'+2
loc extra3=`=r(max)'+3
loc extra4=`=r(max)'+4
set obs `extra4'

qui sum obs
replace significant=1 in `extra1'/`extra2'
replace significant=0 in `extra3'/`extra4'
replace obs=`=r(min)' in `extra1'/`extra4'
replace obs=`=r(max)' in `extra2'/`extra3'

qui sum me_est,detail												/* Marjinal etki tahminlerine dair ölçek değerleri belirlenmektedir. */
loc locut=`r(min)'+(`r(max)'-`r(min)')*1/5 
loc lmedcut=`r(min)'+(`r(max)'-`r(min)')*2/5
loc hmedcut=`r(min)'+(`r(max)'-`r(min)')*3/5
loc hicut=`r(min)'+(`r(max)'-`r(min)')*4/5
loc minest =`r(min)'
loc maxest =`r(max)'

loc colr=`""0 109 44" "65 171 93" "161 217 155" "229 245 224" "242 242 242""' 				/* Renk skalası oluşturulmaktadır. */

twoway (contour me_est x1 x2 if me_est!=., ccuts(`locut' `lmedcut' `hmedcut' `hicut') ccolors(`colr')) ///
(scatter x1 x2 [fw=obs] if significant==0, msymbol(oh) mcolor(black%75) mlwidth(vvthin) msize(vsmall)) /// 
(scatter x1 x2 [fw=obs] if significant==1, msymbol(o) mcolor(black%75) mlwidth(none) msize(vsmall)), ///
xtitle("X{subscript:2}", height(4)) ytitle("X{subscript:1}", height(4)) ztitle("") zlabel(`minest' `locut' `lmedcut' `hmedcut' `hicut' `maxest') ///
legend(off) clegend(title("") ring(1) width(5) height(25)) ///
title("Şekil 5a. Ort. Marjinal Etki Tahminleri (Bileşik Dağılım)") name(Şekil5a, replace)	/* Kontur grafiği çizilmektedir. */

use `DEML', clear
logit

qui sum x2 															/* Koşullandırıcı değişken özetlenmektedir. */
loc mn=r(min)														/* Asgari ve azami değerleri belirlenmekte, ardından grafiklenmek istenen sayıda (10) eşit boyutta alt örnekleme bölünerek bunların ortanca noktaları belirlenmektedir. */
loc mx=r(max)
xtile group_id =x2, nq(10)
egen midpoint=median(x2),by(group_id)

collapse (count) wt=y3, by(midpoint x1 x2 x3)						/* Gözlem sayıları, ortanca noktalar ve model eşitliğinde yer alan tüm kurucu ve bağımsız değişkenlerden mütevellit bir veriseti oluşturulmakta, */
gen x1x2=x1*x2 														/* ardından çarpımlı etkileşim oluşturularak, bu değerler X matrisine aktarılmaktadır. */
putmata wt=wt group_id=midpoint X=(x1 x2 x1x2 x3 1), replace
preserve 
replace x1 =x1+1
replace x1x2=x1*x2
putmata X1=(x1 x2 x1x2 x3 1), replace
restore

mata: dame=me_wt(coef, X, X1, group_id, wt)							/* Dağılım-ağırlıklı marjinal etki tahminleri, yukarıda tanımlanan program yardımıyla hesaplanmaktadır. */

collapse (mean) x1 (median) x3 [fw=wt] 								/* Karşılaştırmak için simülasyon yöntemiyle x1'in ortalama marjinal etkisi, x2'nin 21 farklı değeri için ([0,1] arasında ve 0.05 birimlik artışlarla) yukarıda tanımlanan ilk program yardımıyla hesaplanmaktadır. */
expand 21
gen x2=`mn' + (_n-1)*(`mx'-`mn')/20
gen x1x2=x1*x2

putmata X=(x1 x2 x1x2 x3 1), replace
preserve 
replace x1 =x1+1
replace x1x2=x1*x2
putmata X1=(x1 x2 x1x2 x3 1), replace
restore
 
mata: mem=me_byrow(coef, X, X1)
getmata (mem lbm ubm)=mem 											/* Bu programdan ortalama marjinal etkiler ve bunların güven aralıkları da elde edilebilmektedir. */
getmata (midpoint obs dame_est lb ub)=dame, force

twoway (line mem x2, lpattern(solid) lwidth(thin) lcolor(black%50)) ///
(rline lbm ubm x2, lpattern(shortdash) lwidth(thin) lcolor(black%50)) ///
(scatter dame_est midpoint [fw=obs], msymbol(o) mcolor(maroon%75) mlwidth(none) msize(small)) /// 
(rspike lb ub midpoint, lcolor(maroon%50) lwidth(thin)), ///
yline(0, lcolor(red)) ysca(alt) ylab(-.45(.15).15) legend(order(3 "X{subscript:1}'in Dağılım-Ağırlıklı Ort. Marjinal Etkisi" 1 "Pr(Y{subscript:3}=1 | X{subscript:1}=D{subscript:ort.+1}) - Pr(Y{subscript:3}=1 | X{subscript:1}=D{subscript:ort.})") rows(2) ring(0) pos(7) size(vsmall)) ///
ytitle("X{subscript:1}'in Dağılım-Ağırlıklı Ortalama Marjinal Etkisi") xtitle("X{subscript:2}") title("Şekil 5b. Dağılım-Ağırlıklı Ort. Marjinal Etki Tahminleri") graphregion(margin(l=22)) name(Şekil5b, replace)

gr combine Şekil5a Şekil5b, xcommon 
gr_edit .plotregion1.graph1.yaxis1.title.DragBy 0 -1.25
gr export "Şekil5.jpg", as(jpg) quality(100) replace
